clear
clc
set(0,'DefaultFigureWindowStyle','docked')

%% Notes: The method used is described in Elsby, Michaels, Solon (2009 AEJ, Macroeconomics)

%% read data
% note: all data is seasonally adjusted

data_num = xlsread('match_eff_input_data.xlsx','Data');

yyyymm = data_num(:,1);
L = data_num(:,2);
U = data_num(:,3);
U_less5wks = data_num(:,4);

clear data_num

% format:
% 1. date: yyyymm
% 2. CLF16OV: Civilian Labor Force, Thousands of Persons, Monthly, Seasonally Adjusted
% 3. UNEMPLOY: Unemployment Level, Thousands of Persons, Monthly, Seasonally Adjusted
% 4. UEMPLT5: Number of Civilians Unemployed for Less Than 5 Weeks, Thousands of Persons, Monthly, Seasonally Adjusted

% Feb. 1994: Discontinuity in how CPS measures short-term unemployment 
T_break = 199402; 
Ellsby_multiplier = 1.1549; % See bottom of page 94
U_less5wks(yyyymm>=T_break) = Ellsby_multiplier*U_less5wks(yyyymm>=T_break);

%% monthly job finding rates

T = numel(yyyymm);

F = NaN*zeros(T,1); % job finding probabilities
f = NaN*zeros(T,1); % job finding hazard rate
s = NaN*zeros(T,1); % separation hazard rates

for t = 1:T-1
    F(t) = 1 - ((U(t+1) - U_less5wks(t+1))/U(t));
    f(t) = -log(1-F(t));
end

s0 = 0.03;
options = optimoptions('fsolve','Display','off');
err_s = zeros(T,1);
exitflag_s = zeros(T,1);

for t = 1:T-1
    U_ss = @(s) s*L(t)./(s+f(t));
    seprate_eq = @(s) ( -U(t+1) + U_ss(s) -  ( U(t) -U_ss(s))*exp(- (s +  f(t) )) );
    [s(t), err_s(t), exitflag_s(t)] = fsolve(seprate_eq, s0, options);
end

tolerance = 1e-10;
if max(abs(err_s)) > tolerance
    error('Error in solving for separation rates.');
end

%% plot data

M = mod(yyyymm,100);
Y = (yyyymm - M)/100;
datenums = datenum(Y,M,1);

figure

subplot(1,3,1)
hold on
semilogy(datenums, L)
datetick('x')
title('Civilian labor force')
axis square
recessionplot;

subplot(1,3,2)
hold on
semilogy(datenums, U)
datetick('x')
title('Unemployment level')
axis square
recessionplot;

subplot(1,3,3)
hold on
semilogy(datenums, U_less5wks)
datetick('x')
title('Unemployed for less than 5 weeks')
axis square
recessionplot;

%% plot flow rates

figure
plot(datenums, F)
datetick('x')
title('F')
axis square
recessionplot;

figure
plot(datenums(3:end), (F(1:end-2) + F(2:end-1) + F(3:end))/3)
datetick('x')
title('F, 3 mo ma')
axis square
recessionplot;

figure
plot(datenums, f)
datetick('x')
title('f')
axis square
recessionplot;

figure
plot(datenums(3:end), (f(1:end-2) + f(2:end-1) + f(3:end))/3)
datetick('x')
title('f, 3 mo ma')
axis square
recessionplot;

figure
plot(datenums, s)
datetick('x')
title('s')
axis square
recessionplot;

figure
plot(datenums(3:end), (s(1:end-2) + s(2:end-1) + s(3:end))/3)
datetick('x')
title('s, 3 mo ma')
axis square
recessionplot;

%% save

xlswrite('unemp_flow_rates.xlsx',{'yyyymm','F','f','s'},'Sheet1','A1:D1')
xlswrite('unemp_flow_rates.xlsx',[yyyymm(1:T-1),F(1:T-1),f(1:T-1),s(1:T-1)],'Sheet1',['A2:D',num2str(T)])